arXiv:1501.06434vl [astro-ph.HE] 26 Jan 2015 


Galactic Propagation of Cosmic Rays from Individual 
Supernova Remnants 


Nils Nierstenhoefer, Philipp Graeser, Florian Schuppan and Julia 
Becker Tjus 

Ruhr Universitat Bochum, Universitatsstrafie 150, D-44780 Bochum 
E-mail: nils@tp4.rub.de 


Abstract. It is widely believed that supernova remnants are the best candidate sources for 
the observed cosmic ray flux up to the knee, i.e. up to ~PeV energies. Indeed, the gamma-ray 
spectra of some supernova remnants can be well explained by assuming the decay of neutral 
pions which are created in hadronic interactions. Therefore, fitting the corresponding gamma 
spectra allows us to derive the spectra of cosmic rays at the source which are locally injected 
into our Galaxy. Using these spectra as a starting point, we propagate the cosmic rays through 
the Galaxy using the publicly available GALPROP code. Here, we will present first results on 
the contribution of those SNRs to the total cosmic ray flux and discuss implications. 


1. Introduction 

Many popular theories for the origin of Galactic cosmic rays (CRs) up to the knee-region (~PeV) 
are based on the idea of shock acceleration of charged particles in supernova remnants (SNRs) [1- 
3]. That SNRs indeed provide the required energy budget to be the birthplace of the Galactic 
cosmic ray population has firstly been pointed out in [4]. 

Recent indirect observational techniques as used by the instruments H.E.S.S., MAGIC, 
VERITAS and MILAGRO allow for the observations of very high energy (VHE) gamma-rays 
from SNRs up to ~100 TeV. Including radio and X-ray observations, it is possible to test 
and distinguish models of leptonic or hadronic origin of those gamma-rays as has been done 
systematically for 24 well-identihed SNRs in [5]. In particular, this kind of analysis provides a 
parametrization jp{T) of a potential population of CR protons inside the SNRs. This indirect 
approach to derive jp{T) is useful as the CRs trajectories, masses and energies are altered during 
their Galactic propagation - thus crucial direct information is lost. This is different for the VHE 
gamma-rays which propagate rectilinear with no significant energy losses through the Galaxy. 
For the same reason, the effects of Galactic CR propagation have to be taken into account if one 
aims on testing whether the source spectra jp{T) of SNRs can explain the CR flux as measured 
on Earth. To do so, an altered version of the publicly available GALPROP framework [6] for 
Galactic propagation is used, in which CR nuclei are injected according to the source spectra 
derived in [5]. It needs to be taken into account that these spectra represent currently observed 
SNRs, and thus only a short temporal phase of cosmic ray injection: while a single SNR is 
only active for about 10^ years, cosmic rays diffuse in the Galaxy for about 10® — 10^ years. In 
order to have a time-averaged sample of SNRs, we assume that the CR spectra as derived from 
the currently observed SNRs are a representative subsample of SNR spectra and luminosities. 



We then perform multiple simulations with SNR locations randomly drawn from the spatial 
distribution of [7], which represents a model for the distribution of SNRs in the Galaxies, with 
an enhanced density near the Galactic Genter. 

The main goal of this paper is to describe our current simulation approach and to present 
an exemplary simulation which is focusing on the investigation of the shape of the predicted 
CR proton flux at Earth. In particular, we emphasize the fact that we do not perform any 
optimization of model parameters here and that we present the method of our approach with only 
preliminary physics results in these proceedings. All details of the physics questions, including 
a parameter study will be discussed elsewhere in the future. One of the primary, final goals of 
developing this method is the evaluation of the overall flux normalization. 

Finally, in these proceedings, we present improvements of the simulation approach and their 
impact on the interpretation of the results are outlined. 


2. Source Spectra 

The modeling of the electromagnetic spectrum of 24 SNRs has been performed in [5]. Out of 
those, a significant contribution of gamma-rays from hadronic interactions with a subsequent 
TTo-decay are found for 21 SNRs (see table 3 in [5]). These 21 source spectra jp{T) are the input 
for the analysis which is presented in what follows. 

More precisely, the authors of [5] apply a one zone model for the SNR and provide 
a parametrization of a simple power law for the proton spectra in momentum p = 
— {T + mc^Y, which translates into an energy spectrum of 


jp{T) = Up 
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Here, T, a and m denote the kinetic energy, the normalization constant and the proton mass, 
respectively. A smooth low energy cut off at Tmin = 10 MeV and its high energy counterpart at 
Tmax is included through the factors four and five in equation (1). In this paper, the focus is on 
the CR energy range from ~GeV to ~PeV. Hence, the low energy cut off is not applied. 

It should be noted that for example 274 SNRs are listed in the well known Greens catalogue [8] 
while only ~20 SNRs have been detected in gamma-rays. Correspondingly, the prediction of a 
too low CR-flux is expected here. 

Only SNRs that are observed in MeV-GeV or GeV-TeV are used in the analysis of [5] and 
are included in our analysis. Therefore, the way of selecting sources for the observation with 
e.g. Cherenkov telescopes, may induce a selection bias (e.g. favouring flat, luminous SNRs). 
Therefore, the set of 21 SNRs from [5] is not statistically complete and, consequently, the 
presented results need to be interpreted with adequate care. Future data taken with the 
instruments CTA and HAWC will provide a more complete picture of gamma-ray emitting 
SNRs in the Galaxy and this study can be updated accordingly. 


3. Propagation 

The effects of Galactic propagation on the source spectra jp{T) was simulated using version 54 
of the open source C-|—I—framework GALPROP [6,9,10]. The code was extended to allow for 
particle injection according to jp(T) as parametrized in equation (1). Additionally, a procedure 
was implemented that allows to randomly select a position for a set of SNRs from a predefined 
distribution of matter within the Galaxy. 


3.1. Normalization 

First of all it should be noted that the CR normalization in standard GALPROP applications 
is usually fixed by globally scaling the calculated cosmic ray density such that the observed CR 



flux at Earth is met. In the approach presented in these proceedings, the idea is to fix the 
normalization by the rate of particle injection of the individual SNRs into the Galaxy^. Hence, 
from the technical viewpoint a key ingredient was to rewrite the source spectrum jp(T) in terms 
of the internal GALPROP units. A pragmatic approach to determine the needed conversion 
factor is to compare the calculation of the luminosity L in GALPROP with the corresponding 
integral expression on the basis of jp(T) 

rle9 MeV 

L = aRl^j,c (ITf^TMT). ( 2 ) 

JlO MeV 

Here, Rsnr is the radius of the SNRs. Note that a = 1 corresponds to the case where the CR 
particles stream out of the SNR in radial direction. In contrast to that a =1/2 would arises 
from an averaging assuming that the GRs move in random directions inside the SNRs. The 
aforementioned comparison leads to the following relation between the initial source function 
qi{p{T)) as implemented in GALPROP^ and the spectrum parametrization jp{T) 

iMT)) = Jp(r)- (3) 

Alternatively, the luminosity can be derived from the total energy Etot of protons in the SNR 
via L = Etot It with a time-scale r, representing the distribution of the total energy over the 
total lifetime of the remnant. 


1 Ajr fle9 MeV 

L=-TrRl^K dTTMT). (4) 
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Using this expression for the luminosity, one finds an alternative expressions^ for equation (3) 

( 5 ) 

oLgridT 

Assuming that Etot is the energy of the SNR converted into protons and r is the life time of the 
SNR, L = Etot It can be interpreted as the average luminosity in GRs. 

It is a typical approach in gamma-ray astronomy to derive the total cosmic ray energy budget 
in order to estimate the SNRs possible contribution to the total cosmic ray budget. A back-of 
the envelope calculation predicts that the observed cosmic ray luminosity of the Galaxy (about 
3T0^° erg/s) can be reproduced if on average, 10^*^ erg are going into a single SNR at a Supernova 
rate of 1/(100 years). Here, it is assumed that cosmic rays are injected into the interaction region 
continuously at a constant rate during the lifetime r with the total energy going into cosmic 
rays conserved over time. It is clear that this is a simplifying assumption, as it is known that 
at least the energy spectrum is changing with time, in particular concerning the reduction of 
maximum energy, see e.g. [11]. For energy spectra steeper than the total energy budget is 
dominated by the lower integration threshold, so effects from this temporal development should 
be relatively small. It is also specified in [11] that the total energy of the SNR is decreasing 
with time due to cooling effects. This would mean that, if we assume a constant fraction of the 

^ In the exemplary simulation discussed later in this paper the predicted CR flux is normalized to experimental 
data. The more extensive subject of modeling the CR normalization as a direct consequence of the source spectra 
jp{T) will be discussed in a future paper. 

^ Further details can be found in the GALPROP documentation [9,10] or by directly checking the source code 
file crJuminosity.ee. 

^ Equation (3) has been applied in this paper. 



SNR energy going into cosmic rays at a given time, the actual average luminosity of cosmic rays 
would be underestimated in particular for old remnants. It is not clear, however, if the fraction 
of energy going into cosmic rays is actually constant over time or if it actually decreases with the 
available energy budget. Thus, we judge that in hrst order approximation, it seems reasonable 
to estimate the total energy of the remnant from the given value, keeping in mind the above 
discussed uncertainties. 

In these proceedings, we refrain from discussing the global normalization and normalize the final 
spectrum to the observed flux. The normalizations of the individual SNRs relative to each other 
are fixed using equation 3. 

3.2. Including CR Nuclei 

Although the source spectra taken from [5] are only provided for CR protons, it is possible to 
include CR nuclei in GALPROP simulations (see chapter 5.5 in [9,10]). To do so, the initial 
source function of nuclei ?a(pa) with mass number A and momentum pA is related to qi{pi) by 
the relative abundance X according to 

^ ( 6 ) 

In this context two remarks have to be made: 

(i) Due to the high energy cut off in equation (1) X is not independent of the energy in what 
follows. 

(ii) Including CR nuclei injection, the total energy in hadrons of the SNRs is artihcially 
increased. Note, one can approximately correct the energy budget derived in [5] by down- 
scaling the proton normalization Op in equation (1) appropriately. 

Here, simulations are performed for all nuclei, but the resulting energy spectra are only discussed 
for protons. In the future, the heavy nuclei spectra will be discussed as well in order to investigate 
other questions connected to cosmic ray observations, sources and transport. 

3.3. Random Positions, Mapping to Closest Gridpoint and Averaging 

The idea behind GALPROP is to numerically solve the transport equations that governs the 
propagation of CRs in the Galaxy [6]. This requires the definition of a spatial simulation grid 
on which the particle densities are calculated. Here, we select random position for the SNRs 
and map them to the closest grid point. As the grid size is larger than the extensions of the 
SNRs, this implies that currently the SNRs are treated as point like CR sources. 

We hnally calculate and present average values - e.g. the average proton flux < dF/dT > or the 
Boron to Carbon ratio < B jC > - from N sets of 21 SNRs. It is this averaging which allows 
one to derive a realistic prediction for CR observables. 

Note that the variance of CR observables e.g. Xai{dF/dT) or related quantities cannot be 
conclusively calculated within our current simulation procedure. The reasons for that and and 
a method to overcome this issue in the future is outlined in section 5. 

3 . 4 . GALPROP Settings 

As a starting point we use the standard settings as suggested by the GALPROP web-run 
application web page [9,10]. One change is applied to the standard settings: The size of 
the galaxy has been reduced in our example simulation to -10 kpc< x,y <10 kpc - to allow 
faster simulations. In particular we do not aim for optimizing free model parameters to improve 
the agreement with experimental data at this point. Here, we mainly intend to outline our 
computational approach and to show exemplary results - further results (including more detailed 
parameters studies) will be reported elsewhere in the future. 




Figure 1. 20 out of 100 proton spectra dF/dT of random-position sets of 21 SNRs are 

exemplarily presented in grey. Mean values < dF/dT > (red) and their errors (yellow band) 
are displayed. Blue dots show measured data by CREAM [12], the magenta curve presents the 
observations by PAMELA [13]. Note that that < dF/dT > is normalized to fit the CREAM 
data point at lowest energy. 


4. Exemplary Simulation 

In this section an exemplary simulation is presented. The basic parameters for the underlying 
GALPROP runs are mainly the ones suggested by the GALPROP web-run web page [9,10] (cf. 
section 3). In this example A^=100 sets of 21 SNRs with random positions according to the 
SNR-distribution of [7] have been simulated. In what follows, the bracket-notation (< ... >) 
indicates that the predicted observables that can be compared to the observations are calculated 
as an average of A^=100 simulated sets of SNRs with random positions (cf. section 3 for details). 

The predicted fiux < dF/dT > of cosmic ray protons is compared to experimental data from 
GREAM [12] and PAMELA [13] in Eig. 1. Here, the statistical and systematic experimental 
uncertainties are displayed. Additionally, < dF/dT > has been normalized to CREAM data 
at 3.25 TeV to allow for a comparison of the spectral behaviour. Note that the proper solar 
modulation potential of 550 MV linked with the PAMELA measurements [13] is included in 
Eig. 1 according to [14]. 

We find that the spectral behavior of the measured proton spectrum is well described in the 
energy range of the CREAM experiment, but a discrepancy with the PAMELA data which is 
most obvious at 10 GeV is apparent. As we derive the individual proton spectra from gamma- 
ray data, we expext that best agreement should be present in the statistically well-measured 
range, which is around 1 — 100 GeV for Eermi and 100 GeV-10 TeV for lACTs. Assuming that 
about 10% of the cosmic ray energy is going into photons, this would mean a good description 
between 10 GeV and 100 TeV. The disagreement is therefore not arising from any extrapolation 





























Figure 2. 20 out of 100 -B/C-ratios of random-position sets of 21 SNRs are exemplarily 

presented in grey. Mean values < BjC > (red) and their errors (yellow band) are displayed. 
Blue dots show measured data from AMS-01 [15]. 


effect. 

As a cross-check, we also include the simulation of heavy nuclei and find that the simulated 
< BjC >-ratio is in good agreement with data measured by AMS-01 as can be seen in Fig. 2. 

5. Discussion 

In general the good agreement between < BjC > and experimental data indicates that the used 
simulation setup is well suited to model Galactic CR propagation. Taking into account that 
the observed proton spectrum is well reproduced for energies larger than ~1 TeV, supports the 
famous theory that SNRs are the birthplace of the Galactic GRs up to ~PeV energies. 

These first results show that there is a possible discrepancy between < dF/dT > and the 
observations by PAMELA which is most apparent at kinetic energies of T ~ 10 GeV, using 
GALPROP standard parameters. In the following, we will discuss both statistical and physics 
reasons for the disagreement. 

Statistical reasons may contribute to the disagreement: our current simulation approach 
predicts the average GR observable, e.g. < dF/dT >, but not the proper corresponding variance, 
e.g. Var{dF/dT) (see the discussion later on). However, the variance or alternative statistical 
measures are needed to fully quantify whether the aforementioned discrepancy of < dF/dT > 
is statistically significant. To quantify the agreement between the predicted GR observables 
statistically, a measure of the variance is needed. The variance of the GR observables could be 
calculated and presented here along with the mean value. However, it should be noted that this 
variance has a limited statistical interpretation in the current simulation approach. It merely 
measures the spread of the simulated CR observable which is induced by randomly selecting 


























positions for the 21 SNR in the Galaxy. In particular, this variance would decrease if the 
number of SNRs would be increased. As more than 21 SNRs are expected to contribute to the 
CR flux the variance calculate here would presumably be an upper limit for the true variance. 
In addition, the calculation of the true variance would require to include temporal aspects of 
the SNRs - such as their production rate and lifetime - into our simulation approach. This 
could be done as a generalization of our procedure to include spatial aspects via the selection 
random positions and mapping of SNRs to the closest grid point on the spatial simulation grid 
(cf. subsection 3.3). As similar issues are already addressed in the GALPROP manual [9,10] 
(see chapter 6), an implementation including the temporal aspects of SNRs should be feasible. 
A nearby source contributing to the measured GR flux is another possible explanation for the 
disagreement in the proton spectra. Extending the method to provide a meaningful measure of 
the variance induced by spatial and temporal aspects of the SNRs will help to estimated the 
theoretical probability that such a constellation is realized in nature. 

Possible physics reasons which can contribute to this potential disagreement are the following: 

(i) By using the standard GALPROP settings, we use a diffusion scalar with Kolmogorov-like 
diffusion, leading to a general steepening of the energy spectrum of a distant source by a 
power of ~ 1/3. Measuerments like the one of the Boron-to-Carbon ratio show that the 
scalar approximation would even allow for a larger steepening [16]. This would presumably 
give a better fit at low energies, while a discrepancy would show towards higher energies 
instead. 

(ii) The latter behavior could be explained by missing, flat SNRs, which are not seen yet 
due to missing sensitivity at the highest energies: The SNR sample extracted from [5] is 
not statistically complete (compare section 2). It is expected to have a more complete 
sample of SNRs in the future, with available GTA and HAWG data. Also, future results 
from IceGube, KM3NeT and IceCubeGen2 will presumably help to clarify which part of 
the gamma-ray signal actually is of hadronic nature. Gorrelation studies between gamma- 
rays and ionization signatures can also contribute to identifying those sources that are 
hadronically dominated, see e.g. [17-19]. 

(iii) Another possibility would be a diffusion scalar with a power law behavior broken at around 
100 GeV, as already discussed in [9]. 

(iv) Finally, the full treatment of the diffusion process with a tensor might change the picutre. 
Here, the difficulty is that the structure of the tensor itself is not very well known. First 
investigations of tensor diffusion in the Galaxy are presented inplemented in the DRAGON 
code ([20], so far optimized to the propagation of electrons [21]) and the PICARD code 
([22], published with first test examples). 

We already stressed that the parameters which define the framework for the Galactic 
propagation have not been optimized in our example application. Especially, using a broken 
or steeper power law for the diffusion coefficient as function of energy may provide a handle to 
reproduce the shape of the experimental proton spectrum at T ~GeV. 

6. Conclusion 

Based on the publicly available GALPROP code for Galactic CR propagation, a procedure is 
introduced to predict CR observables for a given set of hadronic source spectra from SNRs. This 
method is based on calculating the average over multiple GALPROP run using sets of SNR with 
positions randomly distributed according to a given spatial distribution. In these proceedings, 
we illustrate the current procedure in which we use the proton spectra of 21 SNRs inferred 
from the corresponding gamma-ray measurements as presented in [5]. While the < B/C > 
ratio is well reproduced in this simulation we find a discrepancy for the predicted proton flux 
< dF/clT >. As we only present an exemplary study with no parameter study at this point. 



we do not draw any conclusions from these findings yet. The aim of these proceedings is the 
presentation of the method itself. In the future, we will present a full parameter study, including 
the discussion of the normalization of the spectrum as well. The latter aspect will contribute to 
answering the question if those SNRs observed in gamma-rays are representative to explain the 
cosmic ray energy budget as observed at Earth. 

In the future, we particularily plan to extend our procedure to (I) include temporal aspects 
of SNR evolution, (II) include additional SNRs and (III) implemented on an inhomogeneous 
simulation grid. 

On longer terms we intend to perform Monte-Carlo (MC) simulations of the propagation of 
Galactic CRs. As a basis, we suggest to use the publicly available CRPropa MC-framework 
to study the propagation of ultra-high energy cosmic rays in extra galactic environments [23]. 
Especially the redesigned object orientated structure of the upcoming version 3.0 of CRPropa 
seems to allow for an easy extensions for the propagation of galactic CRs [24]. With today’s 
computer technologies, a Monte Carlo treatment of Galactic CRs down to T~I0-100 TeV seems 
possible with reasonable run times. For lower energies it may be sufficient to switch to a diffusive 
approximation to avoid the time intensive numerical solution of the equation of motion in the 
galactic magnetic field. 
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